{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "fda4cd37",
   "metadata": {
    "colab_type": "text",
    "id": "view-in-github"
   },
   "source": [
    "<a href=\"https://colab.research.google.com/github/NeuromatchAcademy/course-content/blob/master/tutorials/W0D4_Calculus/W0D4_Tutorial2.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "08842220",
   "metadata": {},
   "source": [
    "# Neuromatch Academy: Week 0, Day 4, Tutorial 2\n",
    "# Differential Equations\n",
    "\n",
    "__Content creators:__ John S Butler, Arvind Kumar\n",
    "\n",
    "__Content reviewers:__  \n",
    "\n",
    "__Production editors:__ "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7d80998b",
   "metadata": {},
   "source": [
    "---\n",
    "# Tutorial Objectives\n",
    "\n",
    "A great deal of neuroscience can be modelled using differential equations, from gating channels to single neurons to a network of neurons to blood flow, to behaviour. A simple way to think about differential equations is they are equations that describe how something changes. \n",
    "\n",
    "The most famous of these in neuroscience is the Nobel Prize winning Hodgkin Huxley equation, which describes a neuron by modelling the gating of each axon. But we will not start there; we will start a few steps back.\n",
    "\n",
    "Differential Equations are mathematical equations that describe how something like population or a neuron changes over time. The reason why differential equations are so useful is they can generalise a process such that one equation can be used to describe many different outcomes.\n",
    "The general form of a first order differential equation is:\n",
    "\n",
    "\\begin{align*}\n",
    "\\frac{d}{dt}y(t)&=f(t,y(t))\\\\\n",
    "\\end{align*}\n",
    "\n",
    "which can be read as \"the change in a process $y$ over time $t$ is a function $f$ of time $t$ and itself $y$\". This might initially seem like a paradox as you are using a process $y$ you want to know about to describe itself, a bit like the MC Escher drawing of two hands painting [each other](https://en.wikipedia.org/wiki/Drawing_Hands). But that is the beauty of mathematics - this can be solved some of time, and when it cannot be solved exactly we can use numerical methods to estimate the answer (see next tutorial). \n",
    "\n",
    "\n",
    "In this tutorial, we will see how __differential equations are motivated by observations of physical responses.__ We will break down the population differential equation, then the integrate and fire model, which leads nicely into raster plots and frequency-current curves to rate models.\n",
    "\n",
    "**Steps:**\n",
    "- Get an intuitive understanding of a linear population differential equation (humans, not neurons)\n",
    "- Visualize the relationship between the change in population and the population\n",
    "- Breakdown the Leaky Integrate and Fire (LIF) differential equation\n",
    "- Code the exact solution of an LIF for a constant input\n",
    "- Visualize and listen to the response of the LIF for different inputs\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "usdpbbBl8uHE",
   "metadata": {
    "cellView": "form",
    "execution": {
     "iopub.execute_input": "2021-06-02T02:46:22.172836Z",
     "iopub.status.busy": "2021-06-02T02:46:22.172281Z",
     "iopub.status.idle": "2021-06-02T02:46:22.192248Z",
     "shell.execute_reply": "2021-06-02T02:46:22.191678Z"
    }
   },
   "outputs": [],
   "source": [
    "#@markdown Tutorial Slides\n",
    "# you should link the slides for all tutorial videos here (we will store pdfs on osf)\n",
    "\n",
    "from IPython.display import HTML\n",
    "HTML('<iframe src=\"\" frameborder=\"0\" width=\"960\" height=\"569\" allowfullscreen=\"true\" mozallowfullscreen=\"true\" webkitallowfullscreen=\"true\"></iframe>')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a5985801",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-02T02:46:22.195936Z",
     "iopub.status.busy": "2021-06-02T02:46:22.195374Z",
     "iopub.status.idle": "2021-06-02T02:46:22.534546Z",
     "shell.execute_reply": "2021-06-02T02:46:22.533962Z"
    }
   },
   "outputs": [],
   "source": [
    "# Imports\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d40abcd8",
   "metadata": {
    "cellView": "form",
    "execution": {
     "iopub.execute_input": "2021-06-02T02:46:22.541421Z",
     "iopub.status.busy": "2021-06-02T02:46:22.539624Z",
     "iopub.status.idle": "2021-06-02T02:46:22.631768Z",
     "shell.execute_reply": "2021-06-02T02:46:22.631057Z"
    }
   },
   "outputs": [],
   "source": [
    "# @title Figure Settings\n",
    "import IPython.display as ipd\n",
    "from matplotlib import gridspec\n",
    "\n",
    "import ipywidgets as widgets  # interactive display\n",
    "%config InlineBackend.figure_format = 'retina'\n",
    "\n",
    "# use NMA plot style\n",
    "plt.style.use(\"https://raw.githubusercontent.com/NeuromatchAcademy/course-content/master/nma.mplstyle\")\n",
    "my_layout = widgets.Layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6XJo_I6_NbJg",
   "metadata": {
    "cellView": "form",
    "execution": {
     "iopub.execute_input": "2021-06-02T02:46:22.653892Z",
     "iopub.status.busy": "2021-06-02T02:46:22.645037Z",
     "iopub.status.idle": "2021-06-02T02:46:22.657355Z",
     "shell.execute_reply": "2021-06-02T02:46:22.656848Z"
    }
   },
   "outputs": [],
   "source": [
    "# @title Plotting Functions\n",
    "\n",
    "def plot_dPdt(alpha=.3):\n",
    "  \"\"\" Plots change in population over time\n",
    "    Args:\n",
    "      alpha: Birth Rate\n",
    "    Returns:\n",
    "      A figure two panel figure\n",
    "      left panel: change in population as a function of population\n",
    "      right panel: membrane potential as a function of time\n",
    "  \"\"\"\n",
    "\n",
    "  with plt.xkcd():\n",
    "    time=np.arange(0, 10 ,0.01)\n",
    "    fig = plt.figure(figsize=(12,4))\n",
    "    gs = gridspec.GridSpec(1, 2)\n",
    "\n",
    "    ## dpdt as a fucntion of p\n",
    "    plt.subplot(gs[0])\n",
    "    plt.plot(np.exp(alpha*time), alpha*np.exp(alpha*time))\n",
    "    plt.xlabel(r'Population $p(t)$ (millions)')\n",
    "    plt.ylabel(r'$\\frac{d}{dt}p(t)=\\alpha p(t)$')\n",
    "\n",
    "    ## p exact solution\n",
    "    plt.subplot(gs[1])\n",
    "    plt.plot(time, np.exp(alpha*time))\n",
    "    plt.ylabel(r'Population $p(t)$ (millions)')\n",
    "    plt.xlabel('time (years)')\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "def plot_V_no_input(V_reset=-75):\n",
    "  \"\"\"\n",
    "    Args:\n",
    "      V_reset: Reset Potential\n",
    "    Returns:\n",
    "      A figure two panel figure\n",
    "      left panel: change in membrane potential as a function of membrane potential\n",
    "      right panel: membrane potential as a function of time\n",
    "  \"\"\"\n",
    "  E_L=-75\n",
    "  tau_m=10\n",
    "  t=np.arange(0,100,0.01)\n",
    "  V= E_L+(V_reset-E_L)*np.exp(-(t)/tau_m)\n",
    "  V_range=np.arange(-90,0,1)\n",
    "  dVdt=-(V_range-E_L)/tau_m\n",
    "\n",
    "  with plt.xkcd():\n",
    "    time=np.arange(0, 10, 0.01)\n",
    "    fig = plt.figure(figsize=(12, 4))\n",
    "    gs = gridspec.GridSpec(1, 2)\n",
    "\n",
    "    plt.subplot(gs[0])\n",
    "    plt.plot(V_range,dVdt)\n",
    "    plt.hlines(0,min(V_range),max(V_range), colors='black', linestyles='dashed')\n",
    "    plt.vlines(-75, min(dVdt), max(dVdt), colors='black', linestyles='dashed')\n",
    "    plt.plot(V_reset,-(V_reset - E_L)/tau_m, 'o', label=r'$V_{reset}$')\n",
    "    plt.text(-50, 1, 'Positive')\n",
    "    plt.text(-50, -2, 'Negative')\n",
    "    plt.text(E_L - 1, max(dVdt), r'$E_L$')\n",
    "    plt.legend()\n",
    "    plt.xlabel('Membrane Potential V (mV)')\n",
    "    plt.ylabel(r'$\\frac{dV}{dt}=\\frac{-(V(t)-E_L)}{\\tau_m}$')\n",
    "\n",
    "    plt.subplot(gs[1])\n",
    "    plt.plot(t,V)\n",
    "    plt.plot(t[0],V_reset,'o')\n",
    "    plt.ylabel(r'Membrane Potential $V(t)$ (mV)')\n",
    "    plt.xlabel('time (ms)')\n",
    "    plt.ylim([-95, -60])\n",
    "\n",
    "    plt.show()\n",
    "\n",
    "## LIF PLOT\n",
    "def plot_IF(t, V,I,Spike_time):\n",
    "  \"\"\"\n",
    "    Args:\n",
    "      t  : time\n",
    "      V  : membrane Voltage\n",
    "      I  : Input\n",
    "      Spike_time : Spike_times\n",
    "    Returns:\n",
    "      figure with three panels\n",
    "      top panel: Input as a function of time\n",
    "      middle panel: membrane potential as a function of time\n",
    "      bottom panel: Raster plot\n",
    "  \"\"\"\n",
    "\n",
    "  with plt.xkcd():\n",
    "    fig = plt.figure(figsize=(12, 4))\n",
    "    gs = gridspec.GridSpec(3, 1,  height_ratios=[1, 4, 1])\n",
    "\n",
    "    # PLOT OF INPUT\n",
    "    plt.subplot(gs[0])\n",
    "    plt.ylabel(r'$I_e(nA)$')\n",
    "    plt.yticks(rotation=45)\n",
    "    plt.hlines(I,min(t),max(t),'g')\n",
    "    plt.ylim((2, 4))\n",
    "    plt.xlim((-50, 1000))\n",
    "\n",
    "    # PLOT OF ACTIVITY\n",
    "    plt.subplot(gs[1])\n",
    "    plt.plot(t,V)\n",
    "    plt.xlim((-50, 1000))\n",
    "    plt.ylabel(r'$V(t)$(mV)')\n",
    "\n",
    "    # PLOT OF SPIKES\n",
    "    plt.subplot(gs[2])\n",
    "    plt.ylabel(r'Spike')\n",
    "    plt.yticks([])\n",
    "    plt.scatter(Spike_time, 1 * np.ones(len(Spike_time)), color=\"grey\", marker=\".\")\n",
    "    plt.xlim((-50, 1000))\n",
    "    plt.xlabel('time(ms)')\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "##    Plotting the differential Equation\n",
    "def plot_dVdt(I=0):\n",
    "  \"\"\"\n",
    "    Args:\n",
    "      I  : Input Current\n",
    "    Returns:\n",
    "      figure of change in membrane potential as a function of membrane potential\n",
    "  \"\"\"\n",
    "\n",
    "  with plt.xkcd():\n",
    "    E_L = -75\n",
    "    tau_m = 10\n",
    "    V = np.arange(-85, 0, 1)\n",
    "    g_L = 10.\n",
    "    fig = plt.figure(figsize=(6, 4))\n",
    "\n",
    "    plt.plot(V,(-(V-E_L) + I*10) / tau_m)\n",
    "    plt.hlines(0, min(V), max(V), colors='black', linestyles='dashed')\n",
    "    plt.xlabel('V (mV)')\n",
    "    plt.ylabel(r'$\\frac{dV}{dt}$')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d88bf487",
   "metadata": {
    "cellView": "form",
    "execution": {
     "iopub.execute_input": "2021-06-02T02:46:22.664362Z",
     "iopub.status.busy": "2021-06-02T02:46:22.663774Z",
     "iopub.status.idle": "2021-06-02T02:46:22.667665Z",
     "shell.execute_reply": "2021-06-02T02:46:22.667165Z"
    }
   },
   "outputs": [],
   "source": [
    "# @title Helper Functions\n",
    "\n",
    "## EXACT SOLUTION OF LIF\n",
    "def Exact_Integrate_and_Fire(I,t):\n",
    "  \"\"\"\n",
    "    Args:\n",
    "      I  : Input Current\n",
    "      t : time\n",
    "    Returns:\n",
    "      Spike : Spike Count\n",
    "      Spike_time : Spike time\n",
    "      V_exact : Exact membrane potential\n",
    "  \"\"\"\n",
    "\n",
    "  Spike = 0\n",
    "  tau_m = 10\n",
    "  R = 10\n",
    "  t_isi = 0\n",
    "  V_reset = E_L = -75\n",
    "  V_exact = V_reset * np.ones(len(t))\n",
    "  V_th = -50\n",
    "  Spike_time = []\n",
    "\n",
    "  for i in range(0, len(t)):\n",
    "\n",
    "    V_exact[i] = E_L + R*I + (V_reset - E_L - R*I) * np.exp(-(t[i]-t_isi)/tau_m)\n",
    "\n",
    "    # Threshold Reset\n",
    "    if V_exact[i] > V_th:\n",
    "        V_exact[i-1] = 0\n",
    "        V_exact[i] = V_reset\n",
    "        t_isi = t[i]\n",
    "        Spike = Spike+1\n",
    "        Spike_time = np.append(Spike_time, t[i])\n",
    "\n",
    "  return Spike, Spike_time, V_exact"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "AcwNWVdsQgLz",
   "metadata": {},
   "source": [
    "---\n",
    "# Section 1: Population Differential Equation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "DIZfUH4rQNSR",
   "metadata": {
    "cellView": "form",
    "execution": {
     "iopub.execute_input": "2021-06-02T02:46:22.674748Z",
     "iopub.status.busy": "2021-06-02T02:46:22.674155Z",
     "iopub.status.idle": "2021-06-02T02:46:22.733362Z",
     "shell.execute_reply": "2021-06-02T02:46:22.732867Z"
    }
   },
   "outputs": [],
   "source": [
    "#@title Video 1: Differential Equations Introduction\n",
    "from IPython.display import YouTubeVideo\n",
    "video = YouTubeVideo(id=\"nWKgnouNzGY\", width=854, height=480, fs=1)\n",
    "print(\"Video available at https://youtu.be/\" + video.id)\n",
    "video"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "286b4298",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "To get an intuitive feel of a differential equations, we will start with a population differential equation, which models the change in population [1], that is human population not neurons, we will get to neurons later. Mathematically it is written like:\n",
    "\\begin{align*}\n",
    "\\\\\n",
    "\\frac{d}{dt}\\,p(t) &= \\alpha p(t),\\\\\n",
    "\\end{align*}\n",
    "\n",
    "where $p(t)$ is the population of the world and $\\alpha$ is a parameter representing birth rate.\n",
    "\n",
    "Another way of thinking about the models is that the equation\n",
    "\\begin{align*}\n",
    "\\\\\n",
    "\\frac{d}{dt}\\,p(t) &= \\alpha p(t),\\\\\n",
    "\\text{can be written as:}\\\\\n",
    "\\text{\"Change in Population\"} &= \\text{ \"Birth rate times Current population.\"}\n",
    "\\end{align*}\n",
    "\n",
    "The equation is saying something reasonable maybe not the perfect model but a good start."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "LRwtiFPVQ3Pz",
   "metadata": {
    "cellView": "form",
    "execution": {
     "iopub.execute_input": "2021-06-02T02:46:22.738749Z",
     "iopub.status.busy": "2021-06-02T02:46:22.738176Z",
     "iopub.status.idle": "2021-06-02T02:46:22.793801Z",
     "shell.execute_reply": "2021-06-02T02:46:22.794274Z"
    }
   },
   "outputs": [],
   "source": [
    "#@title Video 2: The Population Equation\n",
    "video = YouTubeVideo(id=\"yaQdEayD9Nw\", width=854, height=480, fs=1)\n",
    "print(\"Video available at https://youtu.be/\" + video.id)\n",
    "video"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fbad77e6",
   "metadata": {},
   "source": [
    "## Coding Exercise 1: Plot the Change in Population\n",
    "\n",
    "Let's investigate what the rate of change of population looks like as a function of the current population.\n",
    "\n",
    "Write Python code to calculate the change in population `dpdt` with an `alpha=0.3` as a function of population `p`.  Uncomment the plot line and interpret the plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "00bc0d0e",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-02T02:46:22.800288Z",
     "iopub.status.busy": "2021-06-02T02:46:22.799205Z",
     "iopub.status.idle": "2021-06-02T02:46:22.886065Z",
     "shell.execute_reply": "2021-06-02T02:46:22.885402Z"
    }
   },
   "outputs": [],
   "source": [
    "p = np.arange(0, 100, 0.1)\n",
    "\n",
    "\n",
    "########################################################################\n",
    "## TODO for students: calculate the dp/dt\n",
    "## Complete line of code\n",
    "raise NotImplementedError(\"Calculate the dp/dt\")\n",
    "########################################################################\n",
    "\n",
    "dpdt = ...\n",
    "\n",
    "## Uncomment to plot answer\n",
    "fig = plt.figure(figsize=(6, 4))\n",
    "#plt.plot(p, dpdt)\n",
    "plt.xlabel(r'Population $p(t)$ (millions)')\n",
    "plt.ylabel(r'$\\frac{d}{dt}p(t)=\\alpha p(t)$')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8e37083e",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-02T02:46:22.891802Z",
     "iopub.status.busy": "2021-06-02T02:46:22.891132Z",
     "iopub.status.idle": "2021-06-02T02:46:23.363019Z",
     "shell.execute_reply": "2021-06-02T02:46:23.362376Z"
    }
   },
   "outputs": [],
   "source": [
    "# to_remove solution\n",
    "p = np.arange(0, 100, 0.1)\n",
    "\n",
    "with plt.xkcd():\n",
    "\n",
    "  dpdt = 0.3*p\n",
    "  fig = plt.figure(figsize=(6, 4))\n",
    "  plt.plot(p, dpdt)\n",
    "  plt.xlabel(r'Population $p(t)$ (millions)')\n",
    "  plt.ylabel(r'$\\frac{d}{dt}p(t)=\\alpha p(t)$')\n",
    "  plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ffcd1cdd",
   "metadata": {},
   "source": [
    "## Think! 1: Linear Differential Equations and the Population Equation\n",
    "\n",
    "1. Why is the population differential equation is known as a linear differential equation?\n",
    "\n",
    "2. How does population size affect the rate of change of the population?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "il7f6AWkXk0B",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-02T02:46:23.367636Z",
     "iopub.status.busy": "2021-06-02T02:46:23.366983Z",
     "iopub.status.idle": "2021-06-02T02:46:23.370761Z",
     "shell.execute_reply": "2021-06-02T02:46:23.370254Z"
    }
   },
   "outputs": [],
   "source": [
    "# to_remove explanation\n",
    "\n",
    "\"\"\"\n",
    "  1. The plot of $\\frac{dp}{dt}$ is a line, which is why the differential\n",
    "  equation is known as a linear differential equation.\n",
    "\n",
    "  2. As the population increases, the change of population increases. A\n",
    "  population of 20 has a change of 6 while a population of 100 has a change of\n",
    "  30. This makes sense - the larger the population the larger the change.\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5846bbb5",
   "metadata": {},
   "source": [
    "## Section 1.1: Exact Solution of the Population Equation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "uqcubpZNBIGl",
   "metadata": {},
   "source": [
    "### Section 1.1.1: Initial Condition\n",
    "The linear population differential equation is known as an initial value differential equation because we need an initial population value to solve it, so here we will set our initial population to:\n",
    "\n",
    "\\begin{align*}\n",
    "&p(0)=1.\\\\\n",
    "\\end{align*}\n",
    "\n",
    "Different initial conditions will lead to different answers, but they will not change the differential equation. This is one of the strengths of a differential equation. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "oYBK4NWYBL5t",
   "metadata": {},
   "source": [
    "###  Section 1.1.2: Exact Solution\n",
    "To calculate the exact solution of a differential equation, we must integrate both sides and use some mathematical tricks - the harder the equation the more obscure the trick.\n",
    "\n",
    "Given the linear population equation \n",
    "\\begin{align*}\n",
    "\\frac{d}{dt}\\,p(t) &= \\alpha p(t),\\\\\\\\\n",
    "p(0)=1,\\\\\n",
    "\\end{align*}\n",
    "has the exact solution:\n",
    "\\begin{align*}\n",
    "p(t)&=e^{\\alpha t}.\\\\\n",
    "\\end{align*}\n",
    "\n",
    "The exact solution written in words is: \n",
    "\n",
    "\\begin{align*}\n",
    "\\text{\"Population\"}&=\\text{\"grows/declines exponentially as a function time and birth rate\"}.\\\\\n",
    "\\end{align*}\n",
    "\n",
    "Most differential equations do not have a known exact solution, so in the next tutorial on numerical methods we will show how the solution can be estimated.\n",
    "\n",
    "A small aside: a good deal of progress in mathematics was due to mathematicians writing taunting letters to each other saying they had a trick that could solve something better than everyone else. So do not worry too much about the tricks."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1257fcc0",
   "metadata": {},
   "source": [
    "####  Coding Exercise 1.1.2: Exact Solution of the Population Equation\n",
    "Let's consider the population differential equation with a birth rate $\\alpha=0.3$:\n",
    "\n",
    "\\begin{align*}\n",
    "\\frac{d}{dt}\\,p(t) = 0.3 p(t),\\\\\n",
    "\\text{with the initial condition}\\\\\n",
    "p(0)=1.\\\\\n",
    "\\end{align*}\n",
    "\n",
    "It has an exact solution \n",
    "\\begin{align*}\n",
    "\\\\\n",
    "p(t)=e^{0.3 t}.\n",
    "\\end{align*}\n",
    "\n",
    "Code and plot the exact solution of the population differential equation $p(t)$ for $\\alpha=0.3$ for the time period $t=0$ to $t=10$. \n",
    "\n",
    "Use the function `np.exp()` for the exponential."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ac14c49f",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-02T02:46:23.376805Z",
     "iopub.status.busy": "2021-06-02T02:46:23.375181Z",
     "iopub.status.idle": "2021-06-02T02:46:23.384195Z",
     "shell.execute_reply": "2021-06-02T02:46:23.383489Z"
    }
   },
   "outputs": [],
   "source": [
    "t = np.arange(0, 10, 0.1) # Time from 0 to 10 years in 0.1 steps\n",
    "\n",
    "########################################################################\n",
    "## TODO for students: calculate the exact solution of p\n",
    "# Fill out function and remove\n",
    "raise NotImplementedError(\"Calculate the exact solution of p\")\n",
    "########################################################################\n",
    "\n",
    "p = ...\n",
    "\n",
    "## Uncomment to plot answer\n",
    "fig = plt.figure(figsize=(6, 4))\n",
    "#plt.plot(t, p)\n",
    "plt.ylabel('Population (millions)')\n",
    "plt.xlabel('time (years)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eeb717b8",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-02T02:46:23.390430Z",
     "iopub.status.busy": "2021-06-02T02:46:23.389858Z",
     "iopub.status.idle": "2021-06-02T02:46:23.603197Z",
     "shell.execute_reply": "2021-06-02T02:46:23.602259Z"
    }
   },
   "outputs": [],
   "source": [
    "# to_remove solution\n",
    "t = np.arange(0, 10, 0.1) # Time from 0 to 10 years in 0.1 steps\n",
    "\n",
    "with plt.xkcd():\n",
    "\n",
    "  p = np.exp(0.3 * t)\n",
    "\n",
    "  fig = plt.figure(figsize=(6, 4))\n",
    "  plt.plot(t, p)\n",
    "  plt.ylabel('Population (millions)')\n",
    "  plt.xlabel('time (years)')\n",
    "  plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "792841c9",
   "metadata": {},
   "source": [
    "## Section 1.2: Parameters (Birth Rate $\\alpha$)\n",
    "One of the goals when design a differential equation is to make it generalisable. For example, so that the differential equation will give reasonable solutions for different countries with different birth rates $\\alpha$. \n",
    "\n",
    "\n",
    "##  Interactive Demo 1.2:  Interactive Parameter Change\n",
    "Play with the widget to see the relationship between $\\alpha$ and the population differential equation as a function of population (left-hand side), and the population solution as a function of time(right-hand side). Pay close attention to the transision point from positive to negative."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "490c0be0",
   "metadata": {
    "cellView": "form",
    "execution": {
     "iopub.execute_input": "2021-06-02T02:46:23.626756Z",
     "iopub.status.busy": "2021-06-02T02:46:23.626155Z",
     "iopub.status.idle": "2021-06-02T02:46:24.132827Z",
     "shell.execute_reply": "2021-06-02T02:46:24.132177Z"
    }
   },
   "outputs": [],
   "source": [
    "# @markdown Make sure you execute this cell to enable the widget!\n",
    "my_layout.width = '450px'\n",
    "@widgets.interact(\n",
    "    alpha=widgets.FloatSlider(.3, min=-1., max=1., step=.1, layout=my_layout)\n",
    ")\n",
    "def Pop_widget(alpha):\n",
    "  plot_dPdt(alpha=alpha)\n",
    "  plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "Cb6TNCE9ts1k",
   "metadata": {},
   "source": [
    "### Think! 1.2: How do changing parameters of the Population Equation affect the outcome?\n",
    "\n",
    "1. What happens when $\\alpha < 0$?\n",
    "2. What happens when $\\alpha > 0$?\n",
    "3. What happens when $\\alpha = 0$?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "KDgZUHS5YDQp",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-02T02:46:24.140207Z",
     "iopub.status.busy": "2021-06-02T02:46:24.139174Z",
     "iopub.status.idle": "2021-06-02T02:46:24.142317Z",
     "shell.execute_reply": "2021-06-02T02:46:24.141711Z"
    }
   },
   "outputs": [],
   "source": [
    "#to_remove explanation\n",
    "\n",
    "\"\"\"\n",
    "  1. Negative values of alpha result in an exponential decrease.\n",
    "  2. Positive Values of alpha in an exponential increase.\n",
    "  3. Alpha equal to 0 is a unique point known as an equilibrium point when the\n",
    "  dp/dt=0 and there is no change in population. This is known as a stable point.\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c3a78544",
   "metadata": {},
   "source": [
    "## Secton 1.3: Population Differential Equation Summary \n",
    "\n",
    "The population differential equation is an over-simplification and has some very obvious limitations: \n",
    "1. Population growth is not exponential as there are limited number of resources so the population will level out at some point.\n",
    "2. It does not include any external factors on the populations like weather, predators and preys.\n",
    "\n",
    "These kind of limitations can be addressed by extending the model.\n",
    "\n",
    "\n",
    "While it might not seem that the population equation has direct relevance to neuroscience a similar equation is used to describe the accumulation of evidence for decision making known as the Drift Diffusion Model which you will see in more detail in the Linear System day in Neuromatch (W2D2).\n",
    "\n",
    "\n",
    "Another differential equation that is similar to the population equation is the Leaky Integrate and Fire Model which we will explore below."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ae82611f",
   "metadata": {},
   "source": [
    "---\n",
    "# Section 2: The Leaky Integrate and Fire Model\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "hGBNAqnVVY3E",
   "metadata": {
    "cellView": "form",
    "execution": {
     "iopub.execute_input": "2021-06-02T02:46:24.147670Z",
     "iopub.status.busy": "2021-06-02T02:46:24.147066Z",
     "iopub.status.idle": "2021-06-02T02:46:24.204079Z",
     "shell.execute_reply": "2021-06-02T02:46:24.204584Z"
    }
   },
   "outputs": [],
   "source": [
    "#@title Video 3: Leaky Integrate and Fire Model\n",
    "video = YouTubeVideo(id=\"htYJ84aVK3s\", width=854, height=480, fs=1)\n",
    "print(\"Video available at https://youtu.be/\" + video.id)\n",
    "video"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4qgkjHOtWcEL",
   "metadata": {},
   "source": [
    "The Leaky Integrate and Fire Model is a linear differential equation that describes the membrane potential ($V$) of a single neuron which was proposed by Louis Édouard Lapicque in 1907 [2].\n",
    "\n",
    "The subthreshold membrane potential dynamics of a LIF neuron is described by\n",
    "\\begin{align}\n",
    "\\tau_m\\frac{dV}{dt} = -(V-E_L) + R_mI\\,\n",
    "\\end{align}\n",
    "\n",
    "\n",
    "where $\\tau_m$ is the time constant, $V$ is the membrane potential,  $E_L$ is the resting potential, $R_m$ is membrance resistance, and $I$ is the external input current. \n",
    "\n",
    "In the next few sections, we will break down the equation and then build it back up to get an intuitive feel of the different facets of the differential equation."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c05087c6",
   "metadata": {},
   "source": [
    "## Section 2.1: LIF without Input \n",
    "\n",
    "First, we will simplify the equation by removing the input, which gives the equation\n",
    "\n",
    "\\begin{align}\n",
    "\\tau_m\\frac{dV}{dt} &= -V+E_L,\\\\\n",
    "\\end{align}\n",
    "\n",
    "which can be written in words as:\n",
    "\n",
    "\\begin{align}\n",
    "\\begin{matrix}\\text{\"Time constant multiplied by the} \\\\ \\text{change in membrane potential\"}\\end{matrix}&=\\begin{matrix}\\text{\"Minus Current} \\\\ \\text{membrane potential\"} \\end{matrix}+\n",
    "\\begin{matrix}\\text{\"resting potential\"}\\end{matrix}.\\\\\n",
    "\\end{align}\n",
    "\n",
    "\n",
    "The equation can be re-arranged to look even more like the population equation:\n",
    "\n",
    "\\begin{align}\n",
    "\\frac{dV}{dt} &= \\frac{-(V-E_L)}{\\tau_m}.\\\\\n",
    "\\end{align}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e00b7c43",
   "metadata": {},
   "source": [
    "### Coding Exercise 2.1: LIF without Input\n",
    "Code the right hand side `dVdt` of the differential equation,\n",
    "\\begin{align}\n",
    "\\frac{dV}{dt} &= \\frac{-(V-E_L)}{\\tau_m},\\\\\n",
    "\\end{align}\n",
    "with the parameters set as:\n",
    "* `E_L = -75`\n",
    "* `V_reset = -50`\n",
    "* `tau_m = 10`\n",
    "\n",
    "and uncommment the code to plot for 'dVdt' as a function of $V$ between -90 and 0.\n",
    "\n",
    "Take note of where the curve of `dVdt` cuts the dashed line indicating `dVdt = 0`. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f4ea3040",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-02T02:46:24.217285Z",
     "iopub.status.busy": "2021-06-02T02:46:24.211079Z",
     "iopub.status.idle": "2021-06-02T02:46:24.220171Z",
     "shell.execute_reply": "2021-06-02T02:46:24.219547Z"
    }
   },
   "outputs": [],
   "source": [
    "# Parameter definition\n",
    "E_L = -75\n",
    "tau_m = 10\n",
    "\n",
    "# Range of Values of V\n",
    "V = np.arange(-90, 0, 1)\n",
    "\n",
    "########################################################################\n",
    "## TODO for students: Complete line of code for the dVdt\n",
    "# Fill out function and remove\n",
    "raise NotImplementedError(\"Complete line of code for the dVdt\")\n",
    "########################################################################\n",
    "\n",
    "dVdt = ...\n",
    "\n",
    "# Uncomment code below to plot\n",
    "fig = plt.figure(figsize=(6, 4))\n",
    "# plt.plot(V, dVdt)\n",
    "plt.vlines(-75, min(dVdt), max(dVdt), colors='black', linestyles='dashed')\n",
    "plt.text(E_L, max(dVdt) + 1, r'$E_L$')\n",
    "\n",
    "plt.hlines(0, min(V), max(V), colors='black', linestyles='dashed')\n",
    "plt.text(-50, 1, 'Positive')\n",
    "plt.text(-50, -2, 'Negative')\n",
    "plt.xlabel(r'$V(t)$ (mV)')\n",
    "plt.ylabel(r'$\\frac{dV}{dt}=\\frac{-(V-E_L)}{\\tau_m}$')\n",
    "plt.ylim(-8, 2)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8dd9906c",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-02T02:46:24.314153Z",
     "iopub.status.busy": "2021-06-02T02:46:24.295627Z",
     "iopub.status.idle": "2021-06-02T02:46:24.590076Z",
     "shell.execute_reply": "2021-06-02T02:46:24.590536Z"
    }
   },
   "outputs": [],
   "source": [
    "# to_remove solution\n",
    "# Parameter definition\n",
    "E_L = -75\n",
    "tau_m = 10\n",
    "\n",
    "# Range of Values of V\n",
    "V = np.arange(-90, 0, 1)\n",
    "dVdt = -(V - E_L) / tau_m\n",
    "\n",
    "with plt.xkcd():\n",
    "\n",
    "  fig = plt.figure(figsize=(6, 4))\n",
    "  plt.plot(V, dVdt)\n",
    "  plt.hlines(0, min(V), max(V), colors='black', linestyles='dashed')\n",
    "  plt.vlines(-75, min(dVdt), max(dVdt), colors='black', linestyles='dashed')\n",
    "\n",
    "  plt.text(-50, 1, 'Positive')\n",
    "  plt.text(-50, -2, 'Negative')\n",
    "  plt.text(E_L, max(dVdt) + 1, r'$E_L$')\n",
    "  plt.xlabel(r'$V(t)$ (mV)')\n",
    "  plt.ylabel(r'$\\frac{dV}{dt}=\\frac{-(V-E_L)}{\\tau_m}$')\n",
    "  plt.ylim(-8, 2)\n",
    "  plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "XRl1cRt5ppJ1",
   "metadata": {},
   "source": [
    "### Think! 2.1: Effect on Membrane Potential $V$ on the LIF Model\n",
    "\n",
    "1. What is the effect on $\\frac{dV}{dt}$ when $V>-75$ mV?\n",
    "2. What is the effect on $\\frac{dV}{dt}$ when $V<-75$ mV\n",
    "3. What is the effect on $\\frac{dV}{dt}$ when $V=-75$ mV?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "h0MIi-MJYdk5",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-02T02:46:24.596359Z",
     "iopub.status.busy": "2021-06-02T02:46:24.595770Z",
     "iopub.status.idle": "2021-06-02T02:46:24.598970Z",
     "shell.execute_reply": "2021-06-02T02:46:24.599453Z"
    }
   },
   "outputs": [],
   "source": [
    "#to_remove explanation\n",
    "\n",
    "\"\"\"\n",
    "  1. For $V>-75$ mV, the derivative is negative.\n",
    "  2. For $V<-75$ mV, the derivative is positive.\n",
    "  3. For $V=-75$ mV, the derivative is equal to $0$ is and a stable point.\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f7bb8882",
   "metadata": {},
   "source": [
    "### Section 2.1.1: Exact Solution of the LIF model without Input\n",
    "\n",
    "Similar to the population equation, we need an initial membrane potential at time $0$ to solve the LIF model. \n",
    "\n",
    "With this equation \n",
    "\\begin{align}\n",
    "\\frac{dV}{dt} &= \\frac{-(V-E_L)}{\\tau_m}\\,\\\\\n",
    "V(0)&=V_{reset},\n",
    "\\end{align}\n",
    "where is $V_{reset}$ is called the reset potential.\n",
    "\n",
    "The LIF model has the exact solution:\n",
    "\\begin{align*}\n",
    "V(t)=&\\ E_L+(V_{reset}-E_L)e^{\\frac{-t}{\\tau_m}}\\\\\n",
    "\\text{ which can be written as: }\\\\\n",
    "\\begin{matrix}\\text{\"Current membrane} \\\\ \\text{potential}\"\\end{matrix}=&\\text{\"Resting potential\"}+\\begin{matrix}\\text{\"Reset potential minus resting potential} \\\\ \\text{times exponential with rate one over time constant.\"}\\end{matrix}\\\\\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "31910ff4",
   "metadata": {},
   "source": [
    "#### Interactive Demo 2.1.1: Initial Condition $V_{reset}$\n",
    "This exercise is to get an intuitive feel of how the different initial conditions $V_{reset}$ impact the differential equation of the LIF and the exact soluation for the equation:\n",
    "\n",
    "\\begin{align}\n",
    "\\frac{dV}{dt} &= \\frac{-(V-E_L)}{\\tau_m}\\,\\\\\n",
    "\\end{align}\n",
    "with the parameters set as:\n",
    "* `E_L = -75,`\n",
    "* `tau_m = 10.`\n",
    "\n",
    "The panel on the left-hand side plot the change in membrane potential $\\frac{dV}{dt}$ is a function of membrane potential $V$ and right-hand side panel plots the exact solution $V$ as a function of time $t,$ the green dot in both panels is the reset potential $V_{reset}$.\n",
    "\n",
    "Pay close attention to when $V_{reset}=E_L=-75$mV."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "03458759",
   "metadata": {
    "cellView": "form",
    "execution": {
     "iopub.execute_input": "2021-06-02T02:46:24.615557Z",
     "iopub.status.busy": "2021-06-02T02:46:24.606571Z",
     "iopub.status.idle": "2021-06-02T02:46:25.184856Z",
     "shell.execute_reply": "2021-06-02T02:46:25.184235Z"
    }
   },
   "outputs": [],
   "source": [
    "#@markdown Make sure you execute this cell to enable the widget!\n",
    "my_layout.width = '450px'\n",
    "@widgets.interact(\n",
    "  V_reset=widgets.FloatSlider(-77., min=-91., max=-61., step=2,\n",
    "                              layout=my_layout)\n",
    ")\n",
    "\n",
    "def V_reset_widget(V_reset):\n",
    "  plot_V_no_input(V_reset)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "IFbVIjvgraCc",
   "metadata": {},
   "source": [
    "#### Think! 2.1.1: Effect of $V_{reset}$ on the solution\n",
    "\n",
    "\n",
    "1. How does the solution look with initial values of $V_{reset} < -75$?\n",
    "2. How does the solution look with initial values of $V_{reset} > -75$?\n",
    "3. How does the solution look with initial values of $V_{reset} = -75$?\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "IQmxKOuRaVGr",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-02T02:46:25.190849Z",
     "iopub.status.busy": "2021-06-02T02:46:25.190263Z",
     "iopub.status.idle": "2021-06-02T02:46:25.195755Z",
     "shell.execute_reply": "2021-06-02T02:46:25.195269Z"
    }
   },
   "outputs": [],
   "source": [
    "#to_remove explanation\n",
    "\n",
    "\"\"\"\n",
    "  1. Initial Values of $V_{reset} < -75$ result in the solution increasing to\n",
    "  -75mV because $\\frac{dV}{dt} > 0$.\n",
    "  2. Initial Values of $V_{reset} > -75$ result in the solution decreasing to\n",
    "  -75mV because $\\frac{dV}{dt}  < 0$.\n",
    "  3. Initial Values of $V_{reset} = -75$ result in a constant $V = -75$ mV\n",
    "  because $\\frac{dV}{dt}  = 0$ (Stable point).\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5dcd105e",
   "metadata": {},
   "source": [
    "## Section 2.2: LIF with Input\n",
    "We will re-introduce the input $I$ and membrane resistance $R_m$ giving the original equation:\n",
    "\n",
    "\\begin{align}\n",
    "\\tau_m\\frac{dV}{dt} = -(V-E_L) + \\color{blue}{R_mI}\\,\n",
    "\\end{align}\n",
    "\n",
    "the input can be other neurons or sensory information."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9862281c",
   "metadata": {},
   "source": [
    "### Interactive Demo 2.2: The Impact of Input\n",
    "The interactive plot below mainpulates $I$ in the differential equation.\n",
    "With increasing input what would you expect to happen to the solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fc6ce7c7",
   "metadata": {
    "cellView": "form",
    "execution": {
     "iopub.execute_input": "2021-06-02T02:46:25.218948Z",
     "iopub.status.busy": "2021-06-02T02:46:25.218328Z",
     "iopub.status.idle": "2021-06-02T02:46:25.441002Z",
     "shell.execute_reply": "2021-06-02T02:46:25.440478Z"
    }
   },
   "outputs": [],
   "source": [
    "# @markdown Make sure you execute this cell to enable the widget!\n",
    "my_layout.width = '450px'\n",
    "@widgets.interact(\n",
    "  I=widgets.FloatSlider(3., min=0., max=20., step=2,\n",
    "                                layout=my_layout)\n",
    ")\n",
    "\n",
    "def Pop_widget(I):\n",
    "  plot_dVdt(I=I)\n",
    "  plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "QDmbE4Usr0k1",
   "metadata": {},
   "source": [
    "### Think! 2.2: Effect of increasing Input\n",
    "\n",
    "1. As $I$ increases, less and less of $\\frac{dV}{dt}$ is below 0. How would this impact the solution?\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "HSNtznIwY6wA",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-02T02:46:25.447603Z",
     "iopub.status.busy": "2021-06-02T02:46:25.446932Z",
     "iopub.status.idle": "2021-06-02T02:46:25.450662Z",
     "shell.execute_reply": "2021-06-02T02:46:25.450050Z"
    }
   },
   "outputs": [],
   "source": [
    "#to_remove explanation\n",
    "\n",
    "\"\"\"\n",
    "\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "968432a1",
   "metadata": {},
   "source": [
    "### Section 2.2.1: LIF Exact Solution\n",
    "\n",
    "The LIF with a constant input is a linear differential equation has a known exact solution:\n",
    "\\begin{align*}\n",
    "V(t)=&\\ E_L+R_mI+(V_{reset}-E_L-R_mI)e^{\\frac{-t}{\\tau_m}}\\\\\n",
    "\\text{which is written as:}\\\\\n",
    "\\begin{matrix}\\text{\"Current membrane} \\\\ \\text{potential\"}\\end{matrix}=&\\text{\"Resting potential\"}+\\begin{matrix}\\text{\"Reset potential minus resting potential} \\\\ \\text{times exponential with rate one over time constant.\" }\\end{matrix}\\\\\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "91b7c87a",
   "metadata": {},
   "source": [
    "#### Coding Exercise 2.2.1: LIF Exact Solution\n",
    "Code the exact solution `V`:\n",
    "$$\n",
    "V(t)= E_L+RI+(V_{reset}-E_L-RI)e^{\\frac{-t}{\\tau_m}},\\\\\n",
    "$$\n",
    "with the parameters set as:\n",
    "* `V_reset = -75,`\n",
    "* `E_L = -75,`\n",
    "* `tau_m = 10,`\n",
    "* `R_m = 10,`\n",
    "* `I = 10,`\n",
    "of the LIF model and uncomment the code to plot the result.\n",
    "\n",
    "Ask yourself, does the result make biological sense?\n",
    "If not, what would you change?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9797df4d",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-02T02:46:25.459899Z",
     "iopub.status.busy": "2021-06-02T02:46:25.455671Z",
     "iopub.status.idle": "2021-06-02T02:46:25.462733Z",
     "shell.execute_reply": "2021-06-02T02:46:25.462133Z"
    }
   },
   "outputs": [],
   "source": [
    "dt = 0.5\n",
    "t_rest = 0\n",
    "\n",
    "t = np.arange(0, 1000, dt)\n",
    "tau_m = 10\n",
    "R_m = 10\n",
    "V_reset = E_L = -75\n",
    "\n",
    "I = 10\n",
    "\n",
    "########################################################################\n",
    "## TODO for students: Complete line of code for the exact Solution V\n",
    "# Fill out function and remove\n",
    "raise NotImplementedError(\"Complete line of code for the exact Solution V\")\n",
    "########################################################################\n",
    "\n",
    "V = ...\n",
    "\n",
    "# Uncomment below to plot the result\n",
    "fig = plt.figure(figsize=(6, 4))\n",
    "#plt.plot(t,V)\n",
    "plt.ylabel('V (mV)')\n",
    "plt.xlabel('time (ms)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "db5816d4",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-02T02:46:25.534718Z",
     "iopub.status.busy": "2021-06-02T02:46:25.532642Z",
     "iopub.status.idle": "2021-06-02T02:46:25.704381Z",
     "shell.execute_reply": "2021-06-02T02:46:25.703488Z"
    }
   },
   "outputs": [],
   "source": [
    "# to_remove solution\n",
    "dt = 0.5\n",
    "t_rest = 0\n",
    "\n",
    "t = np.arange(0, 1000, dt)\n",
    "\n",
    "tau_m = 10\n",
    "R_m = 10\n",
    "V_reset = E_L = -75\n",
    "\n",
    "I = 10\n",
    "\n",
    "V = E_L + R_m*I + (V_reset - E_L - R_m*I) * np.exp(-(t)/tau_m)\n",
    "\n",
    "with plt.xkcd():\n",
    "\n",
    "  fig = plt.figure(figsize=(6, 4))\n",
    "  plt.plot(t,V)\n",
    "  plt.ylabel('V (mV)')\n",
    "  plt.xlabel('time (ms)')\n",
    "  plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0a31b8d1",
   "metadata": {},
   "source": [
    "## Section 2.3: Maths is one thing, but biology matters\n",
    "While the mathematics of the exact solution is exact, it is not biologically valid as neurons spike and defintely do not plateau at a very positive value.\n",
    "\n",
    "To model the firing of a spike, we must have a threshold voltage $V_{th}$ such that if the voltage $V(t)$ goes above it, the neuron spikes\n",
    "$$V(t)>V_{th}.$$\n",
    "We must record the time of spike $t_{isi}$  and count the number of spikes\n",
    "$$t_{isi}=t, $$\n",
    "$$𝑆𝑝𝑖𝑘𝑒=𝑆𝑝𝑖𝑘𝑒+1.$$\n",
    "Then reset the membrane voltage $V(t)$\n",
    "$$V(t_{isi} )=V_{Reset}.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9f12c68e",
   "metadata": {},
   "source": [
    "To take into account the spike the exact solution becomes:\n",
    "\\begin{align*}\n",
    "V(t)=&\\ E_L+R_mI+(V_{reset}-E_L-R_mI)e^{\\frac{-(t-t_{isi})}{\\tau_m}},&\\qquad V(t)<V_{th} \\\\\n",
    "V(t)=&V_{reset},&\\qquad V(t)>V_{th}\\\\\n",
    "Spike=&Spike+1,&\\\\\n",
    "t_{isi}=&t,\\\\\n",
    "\\end{align*}\n",
    "while this does make the neuron spike, it introduces a discontinuity which is not as elegant mathematically as it could be, but it gets results so that is good."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "acc43e8a",
   "metadata": {},
   "source": [
    "### Interactive Demo 2.3.1: Listen to the Spikes\n",
    "This exercise show the relationship between firing rate and the Input for exact solution `V` of the LIF:\n",
    "$$\n",
    "V(t)=\\ E_L+R_mI+(V_{reset}-E_L-R_mI)e^{\\frac{-(t-t_{isi})}{\\tau_m}},\n",
    "$$\n",
    "with the parameters set as:\n",
    "* `V_reset = -75,`\n",
    "* `E_L = -75,`\n",
    "* `tau_m = 10,`\n",
    "* `R_m = 10.`\n",
    "\n",
    "\n",
    "Below is a figure with with three panels; \n",
    "* the top panel is the input, $I,$\n",
    "* the middle panel is the membrane potential $V(t)$ to illustrative the spike $V(t)$ is set to $0$ and then reset to $-75$ mV, \n",
    "* the bottom panel is the raster plot with each dot indicating a spike.\n",
    "\n",
    "As electrophysiologist normally listen to spikes when conducting experiments, so listen to the music of different firing rates by changing the input value $I.$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "55785e26",
   "metadata": {
    "cellView": "form",
    "execution": {
     "iopub.execute_input": "2021-06-02T02:46:25.727624Z",
     "iopub.status.busy": "2021-06-02T02:46:25.726535Z",
     "iopub.status.idle": "2021-06-02T02:46:26.331683Z",
     "shell.execute_reply": "2021-06-02T02:46:26.331155Z"
    }
   },
   "outputs": [],
   "source": [
    "# @markdown Make sure you execute this cell to be able to hear the neuron\n",
    "I = 3\n",
    "t = np.arange(0, 1000, dt)\n",
    "Spike, Spike_time, V = Exact_Integrate_and_Fire(I, t)\n",
    "\n",
    "plot_IF(t, V, I, Spike_time)\n",
    "ipd.Audio(V, rate=len(V))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "54d5af0e",
   "metadata": {},
   "source": [
    "### Interactive Demo 2.3.2: Input on Spikes\n",
    "Manipulate the input into the LIF to see the impact of input on the firing pattern (rate)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ebf27a22",
   "metadata": {
    "cellView": "form",
    "execution": {
     "iopub.execute_input": "2021-06-02T02:46:26.354077Z",
     "iopub.status.busy": "2021-06-02T02:46:26.353234Z",
     "iopub.status.idle": "2021-06-02T02:46:26.859958Z",
     "shell.execute_reply": "2021-06-02T02:46:26.860523Z"
    }
   },
   "outputs": [],
   "source": [
    "# @markdown Make sure you execute this cell to enable the widget!\n",
    "my_layout.width = '450px'\n",
    "@widgets.interact(\n",
    "  I=widgets.FloatSlider(3, min=2.0, max=4., step=.1,\n",
    "                              layout=my_layout)\n",
    ")\n",
    "\n",
    "def Pop_widget(I):\n",
    "  Spike, Spike_time, V = Exact_Integrate_and_Fire(I, t)\n",
    "  plot_IF(t, V, I, Spike_time)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "Rjdocu_Zs9iH",
   "metadata": {},
   "source": [
    "### Think! 2.3: Is the model biologically valid?\n",
    "\n",
    "1. What is the effect of $I$ on spiking?\n",
    "2. Is this biologically valid?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "GjwI3QF3ZLtY",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-02T02:46:26.873557Z",
     "iopub.status.busy": "2021-06-02T02:46:26.872996Z",
     "iopub.status.idle": "2021-06-02T02:46:26.878747Z",
     "shell.execute_reply": "2021-06-02T02:46:26.878226Z"
    }
   },
   "outputs": [],
   "source": [
    "#to_remove explanation\n",
    "\n",
    "\"\"\"\n",
    "  1. As $I$ increases, the number of spikes increases.\n",
    "  2.\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c3cf4d7b",
   "metadata": {},
   "source": [
    "## Section 2.4 Firing Rate as a function of Input\n",
    "\n",
    "The firing frequency of a neuron plotted as a fucntion of current is called an input-output curve (F–I curve). It is also known as a transfer function, which you came accross in the previous tutorial. This function is one of the starting points for the rate model, which extends from modelling single neurons to collections of neurons. \n",
    "\n",
    "By fitting this to a function, we can start to generalise the firing pattern of many neurons, which can be used to build rate models which is discussed later in Neuromatch. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "96bcf5e7",
   "metadata": {
    "cellView": "form",
    "execution": {
     "iopub.execute_input": "2021-06-02T02:46:26.937864Z",
     "iopub.status.busy": "2021-06-02T02:46:26.923285Z",
     "iopub.status.idle": "2021-06-02T02:46:27.258877Z",
     "shell.execute_reply": "2021-06-02T02:46:27.259358Z"
    }
   },
   "outputs": [],
   "source": [
    "# @markdown *Execture this cell to visualize the FI curve*\n",
    "I_range = np.arange(2.0, 4.0, 0.1)\n",
    "Spike_rate = np.ones(len(I_range))\n",
    "\n",
    "for i, I in enumerate(I_range):\n",
    "  Spike_rate[i], _, _ = Exact_Integrate_and_Fire(I, t)\n",
    "\n",
    "with plt.xkcd():\n",
    "  fig = plt.figure(figsize=(6, 4))\n",
    "  plt.plot(I_range,Spike_rate)\n",
    "  plt.xlabel('Input Current (nA)')\n",
    "  plt.ylabel('Spikes per Second (Hz)')\n",
    "  plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "58fdcbe8",
   "metadata": {},
   "source": [
    "---\n",
    "## Summary of LIF model\n",
    "\n",
    "The LIF model is a very nice differential equation to start with in computational neuroscience as it has been used as a building block for many papers that simulate neuronal response.\n",
    "\n",
    "__Strenghts of LIF model:__\n",
    "+ Has an exact solution;\n",
    "+ Easy to interpret;\n",
    "+ Great to build network of neurons.\n",
    "\n",
    "__Weaknesses of the LIF model:__\n",
    "- Spiking is a discontinuity;\n",
    "- Abstraction from biology;\n",
    "- Cannot generate different spiking patterns.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "697364ff",
   "metadata": {},
   "source": [
    "---\n",
    "# Take Home Message"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "_Jmpnq0mSihx",
   "metadata": {
    "cellView": "form",
    "execution": {
     "iopub.execute_input": "2021-06-02T02:46:27.264697Z",
     "iopub.status.busy": "2021-06-02T02:46:27.263974Z",
     "iopub.status.idle": "2021-06-02T02:46:27.332154Z",
     "shell.execute_reply": "2021-06-02T02:46:27.332630Z"
    }
   },
   "outputs": [],
   "source": [
    "#@title Video 4: Take Home Message and Outro\n",
    "video = YouTubeVideo(id=\"TcKWdBD3jcA\", width=854, height=480, fs=1)\n",
    "print(\"Video available at https://youtu.be/\" + video.id)\n",
    "video"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "x0H0zhTNWlop",
   "metadata": {},
   "source": [
    "In this tutorial, we have have seen two differential equations, the population differential equations and the leaky integrate and fire model.\n",
    "\n",
    "\n",
    "We learned about:\n",
    "* The motivation for differential equations.\n",
    "* An intuitive relationship between the solution and the form of the differential equation.\n",
    "* How different parameters of the differential equation impact the solution.\n",
    "* The strengths and limitations of the simple differential equations.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eaec6994",
   "metadata": {},
   "source": [
    "---\n",
    "# Links to Neuromatch Days\n",
    "\n",
    "Differential equations turn up in a number of different Neuromatch days:\n",
    "* The LIF model is discussed in more details in Model Types (Week 1 Day 1) and Real Neurons (Week 2 Day 3).\n",
    "* Drift Diffusion model which is a differential equation for decision making is discussed in Linear Systems (Week 2 Day 2).\n",
    "* Systems of differential equations are discussed in Linear Systems (Week 2 Day 2) and Dynamic Networks (Week 2 Day 4).\n",
    "\n",
    "\n",
    "---\n",
    "# References \n",
    "1. Lotka, A. L, (1920)   Analytical note on certain rhythmic  relations inorganic  systems.Proceedings of the National Academy of Sciences,6(7):410–415,1920.\n",
    "\n",
    "2. Brunel N, van Rossum MC. Lapicque's 1907 paper: from frogs to integrate-and-fire. Biol Cybern. 2007 Dec;97(5-6):337-9. doi: 10.1007/s00422-007-0190-0. Epub 2007 Oct 30. PMID: 17968583.\n",
    "\n",
    "\n",
    "\n",
    "# Bibliography\n",
    "1. Dayan, P., & Abbott, L. F. (2001). Theoretical neuroscience: computational and mathematical modeling of neural systems. Computational Neuroscience Series.\n",
    "2. Strogatz, S. Nonlinear dynamics and chaos: with applications to physics, biology, chemistry, and engineering (studies in nonlinearity), Westview Press; 2 edition (29 July 2014)\n",
    "\n",
    "## Supplemental Popular Reading List\n",
    "1. Lindsay, G. (2021). Models of the Mind: How Physics, Engineering and Mathematics Have Shaped Our Understanding of the Brain. Bloomsbury Publishing.\n",
    "2. Strogatz, S. (2004). Sync: The emerging science of spontaneous order. Penguin UK.\n",
    "\n",
    "## Popular Podcast\n",
    "1. Strogatz, S. (Host). (2020-), Joy of X  https://www.quantamagazine.org/tag/the-joy-of-x/ Quanta Magazine\n"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "include_colab_link": true,
   "name": "W0D4_Tutorial2",
   "provenance": [],
   "toc_visible": true
  },
  "kernel": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
